Quantitative modeling of laser speckle imaging 
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We have analyzed the image formation and dynamic properties in laser speckle imaging (LSI) both experimentally and with 
Monte-Carlo simulation. We show for the case of a liquid inclusion that the spatial resolution and the signal itself are both signif- 
icantly affected by scattering from the turbid environment. Multiple scattering leads to blurring of the dynamic inhomogeneity as 
detected by LSI. The presence of a non-fluctuating component of scattered light results in the significant increase in the measured 
image contrast and complicates the estimation of the relaxation time. We present a refined processing scheme that allows a correct 
estimation of the relaxation time from LSI data. 
© 2008 Optical Society of America 



Laser speckle imaging (LSI) is an efficient and simple 
method for full-field monitoring of dynamics in heteroge- 
neous media. 1 It is widely used in biomedical imaging of 
blood flow 1-5 since it provides access to physiological pro- 
cesses in vivo with excellent temporal and spatial resolution. 

In a more general context LSI can be considered a simpli- 
fied version of the dynamic light scattering (DLS) approach, 
which analyzes the temporal intensity fluctuations of scat- 
tered laser light in order to derive the microscopic proper- 
ties of the scatterers position and motion. In LSI an image 
of dynamic heterogeneities is obtained by analysing the lo- 
cal speckle contrast in the image plane. The local contrast K 
is defined as the ratio of the standard deviation of the inten- 
sity fluctuations to the mean intensity: 1 K = <T area /(I) area - In 
the absence of scatterers motion the contrast takes a value of 
K = 1, while motion blurs the speckle pattern and thus the 
contrast decreases. Therefore the value K can be used to esti- 
mate the time scale of intensity fluctuations and local scatter 
motion. 1 However the quantitative interpretation of the LSI 
data is not straightforward. Multiple scattering of light can 
influence the apparent size of the object due to diffuse blur- 
ring. Access to the local dynamic properties, such as blood 
flow or Brownian motion, is complicated by the complex in- 
terplay between the measured contrast and the full fluctuation 
spectrum of scattered light. The latter is usually not accessible 
if technical simplicity of LSI is to be preserved. 

In this letter we present a quantitative approach to analyze 
LSI images. We address the problem of spatial resolution and 
blurring due to multiple scattering via model experiments and 
Monte-Carlo simulations. Further we will demonstrate, that 
previous attempts to relate quantitatively LSI images to the 
microscopic motion have been hampered by an incorrect data 
analysis. We present a refined processing scheme to access 
this information from a standard LSI experiment. The main 
new element of our analysis is to take into account quantita- 
tively the contribution of the non-fluctuating part of scattered 



intensity. We furthermore suggest a simple experimental pro- 
cedure that allows to access this important quantity. 

Model experiments have been carried out using a homoge- 
neous block of solid Teflon and a home-made heterogeneous 
sample. This medical phantom mimics a liquid inclusion in 
solid tissue. It is obtained by milling a cylindrical hole of di- 
ameter D — 3 mm in a block of solid Teflon. A layer of vari- 
able thickness 0.45 - 2.1 mm separates the cylindrical inclu- 
sion from the interface that is imaged. The void is filled with a 
dispersion of 710 nm polystyrene particles in water. The par- 
ticle concentration is adjusted to match the optical properties 
of the liquid to the solid (volume fraction ca. 1 .3%; scattering 
coefficient fi s = 36 mm -1 , mean free path /* ss 277 pirn and 
anisotropy factor g * 0.9). Thus our sample does not show 
any static scattering differences. The imaging setup has been 
described in detail in Ref. 6. Briefly, the cylindric inclusion 
with flat base of 3 mm diameter oriented to the surface is 
imaged with a CCD camera (PCO Pixelfly, Germany) with 
a standard camera objective lens (/ = 50 mm). The sample 
is illuminated with a diode laser (wavelength 785 nm, max. 
50 mW). The beam is expanded by a slow-rotating ground 
glass in order to reduce statistical noise. 6 For each depth the 
contrast as a function of radial distance r from the inclusion 
center is determined. The results of this procedure for differ- 
ent depths of the inclusion are shown in Fig.^ 

For the Monte-Carlo simulation we followed the photon 
packet approach of light propagation in a turbid medium. 7 
The Henyey-Greenstein phase function was used based on an 
average scattering angle (cos 9). 1 The degree of polarization 
was assumed to decay exponentially. 8 The reflections and re- 
fractions on the boundaries were treated according to Fresnel 
formulas. As the refractive index of the Teflon (n = 1.35) is 
very close to water (« = 1.33) interactions with boundaries 
of dynamic inclusion were not considered. The photon pack- 
ets back-reflected from the sample where registered within a 
numerical aperture 0.17. Only depolarized photons were reg- 
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Fig. 1 . Contrast K as a function of distance r from the center 
of the inclusion for different depths d. Curves represent sim- 
ulation results, symbols - experimental results. Vertical line 
shows the inclusion boundary. Inset: Measured width of the 
liquid solid-boundary as a function of depth (• - simulation 
data, □ - experiment, line is a guide for the eye). 

istered to match the experimental conditions. The field auto- 
correlation function (ACF) gi(r) of the scattered light was 
determined from the photon packet history as explained in 
Ref. 9. 

In order to quantify the effect of image blurring caused by 
multiple scattering we determine the standard deviation of 
the contrast gradient <r gra d . As can be seen in Fig. ^ start- 
ing from the smallest depths the apparent width is increased. 
This increase is of the order of a few P, which is the rele- 
vant length scale for an incident photon to propagate laterally. 
The smearing increases with depth until for large depths the 
width becomes of the order of the depth, as expected for dif- 
fuse light propagation. 10 Our results show the importance of 
diffuse blurring in the image formation already for small and 
moderate depths and can thus provide important guidelines 
for the analysis of actual experiments in biomedical imag- 

1-5 

ing. 

A typical set of simulated correlation functions is presented 
in the Fig|2] If gi(r) is known, then the speckle contrast K of 
the time-integrated speckle fluctuations can be obtained via 
the following relation 11 ' 12 

* 2 = tS- 1 = 1F f \gi(r)\ 2 (l-r/T)dT, (1) 
v) 1 Jo 

where T is the integration time of the detector and f3 the coher- 
ence factor of the detection optics. It is important to note that 
this equation is different to the traditionally used expression 
of Fercher and Briers. 1 As reaffirmed by Durian and cowork- 
ers the correct expression, eq. Q, has to take into account the 
effect of triangular averaging of the correlation function. 1 1 
As shown in Fig. [2 the resulting values of the contrast are 
found in quantitative agreement with the experimental data 
without any adjustable parameter. 

A comparison of Fig.[fland Fig.[5]immediately reveals the 
main limitation of LSI. A single contrast value K is recorded 



Fig. 2. Simulated correlation functions gi(r) for the inclusion 
depth d — 100 ymi. Solid lines show the fits with DWS the- 
ory with baseline. The relaxation times To for the fits were 
obtained by the contrast K inversion (see Fig.[3J. 



in practice as compared to the complex correlation function 
characterizing the full spectrum of intensity fluctuations and 
thus particular care has to be taken which and how informa- 
tion can be extracted. 

In our work we have access to both K and g\ (?) and thus we 
are able to analyze the LSI detection process quantitatively. 
We first analyze the characteristic features of the correlation 
functions presented in Fig. [2] Overall the functions are well 
described by the usual stretched -exponential form derived for 
diffusing-wave spectroscopy (DWS) of a colloid suspension: 
gi(r) = exp(-y V6t/to), where in our case relaxation time 
to = 1/D£q characterizes Brownian motion (diffusion coef- 
ficient D) and y » 2 is a constant. 13, 14 It is worthwhile to 
note that for the case of blood flow the relaxation time to is 
inversely proportional to the root mean square of the flow ve- 
locity . 

The presence of a non-fluctuating static scattering part, 
however, leads to a non-zero baseline. Since the relative 
amount of the static light scattering increases with the dis- 
tance r from the center, the plateau of the ACF increases 
as well. It has been noted previously that a non-zero base- 
line significantly influences or even dominates the resulting 
values of the contrast K. 1 Nevertheless it is rather common 
in the biomedical community to convert the obtained contrast 
values directly to relaxation times neglecting contributions of 
static scattering. 2 ' 4,5, 15-17 

If the detected light is composed of dynamic and static 
components: E(t) = Edit) + E s {t) the measured field corre- 
lation function is defined as following: 18 

Sl(T)=p|£ld(T)| + (l-p), (2) 

where p = (Id) /((Id) + (Is)) characterizes the fluctuating part 
of the detected light intensity and gi&(j) is the field cor- 
relation function of the dynamic component. The resulting 
contrast value can be calculated by substituting eq. @ into 
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Fig. 3. Inverse of the relaxation time 1/tq as a function of 
depth obtained by directly converting the contrast K based 
on Eq[2 neglecting the static part: • - simulation data, A - 
experiment and using the correct procedure based on Eq. |5J 
- simulation and O - experiment. Horizontal dashed line - 
actual value of the relaxation time. 



eq. CJ: 
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/3[p 2 K 2 d + 2p(l-p)K 2 d + (l-p) 2 ], (3) 

= l/Tf Q T \g ld (T)\ 2 (l-T/T)dT and K 2 ld = 

1 /T J \g\ a (t)| (1 - t/T) dr are the normalized variances of 
the intensity and field fluctuations of the dynamic part of the 
detected signal with K^d < Kid- The last term in eq. is 
time independent and reflects the impact of the static scatte- 
ring contribution. Thus the contrast K„, measured in an LSI 
experiment is given by the dynamic contrast of intensity fluc- 
tuations K2d, the dynamic contrast of the field fluctuations 
K\ d and p. 

In an actual LSI experiment an additional processing step 
has to be introduced in order to separate the dynamic and 
static part. The camera exposure time T is usually larger 
compared to the relaxation time tq related to blood flow and 
subsequent frames are separated by a period At larger than 
T. Since At > T > to two sequential frames are free of 
the dynamic component of interest and thus the time and 
space dependent static or pseudo-static component can be 
easily found by cross-correlating sequential frames p (t, r) = 

1 - [{I (t, r)I(t + At, r)>,/<7 (/, r)) 2 -if' 2 . 

Fig. shows a comparison of the correct inversion pro- 
cedure based on eq. as compared to the one neglect- 
ing static scattering. 12 While the correct procedure produces 
values comparable to the relaxation time of Brownian motion 
to the direct conversion can deviate by several orders of mag- 
nitude. 

In conclusion, we studied the image blurring of an object 
buried in turbid media and found that the resolution of the ob- 
tained images can be affected significantly by multiple scat- 
tering. We furthermore introduced a model that reflects the 



impact of the static scattering on the quantitative interpreta- 
tion of LSI images. A simple procedure has been suggested 
to perform a quantitative analysis in actual LSI experiments. 
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